Anisotropic vortex lattice structures in the FeSe superconductor 
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In the recent work by Song et alri-, the scanning tunneling spectroscopy experiment in the stoichio- 
metric FeSe reveals evidence for nodal superconductivity and strong anisotropy. The nodal structure 
can be explained with the extended s-wave pairing structure with the mixture of the s^2^y2 and 
s^2y2 pairing symmetries. We calculate the anisotropic vortex structure by using the self-consistent 
Bogoliubov-de Gennes mean-field theory. In considering the absence of magnetic ordering in the 
FeSe at the ambient pressure, orbital ordering is introduced, which breaks the C4 lattice symmetry 
down to C2, to explain the anisotropy in the vortex tunneling spectra. 

PACS numbers: 74.20.Rp, 71.10.Fd, 74.20.Mn 
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I. INTRODUCTION 

Since the first iron-based layered superconductor 
La(0i_a;F2:)FeAs had been discovered^, the family of 
iron-based superconductors have brought huge impact in 
condensed matter physics community. These novel mate- 
rials exhibit similar phase diagrams compared to high-Tc 
cuprates. The parent compound LaOFeAs has an antifer- 
romagnetic spin-density-wave order—, and, upon doping, 
superconductivity appears. Although correlation effects 
are weaker in iron-based superconductors than those in 
high-Tc cuprates, novel features arise from the multi- 
orbital degree of freedom. The orbital band structures 
play fundamental roles in determining the Fermi surface 
configurations and pairing structures. 

Understanding pairing symmetries is one of the most 
important issues in the study of iron-based supercon- 
ductors. Based on various experimental^^- and theoret- 
ical works the nodeless Sj:2j^2-wave pairing has been 
proposed. In momentum space, the Fermi surfaces of 
many iron-based superconductors consist of hole pock- 
ets around the F-point, and electron pockets around the 
two M-points. The signs of the pairing order param- 
eters on electron and hole Fermi surfaces are opposite. 
The nodal lines of the gap function have no intersections 
with Fermi surfaces, thus the Sa,2y2-pairing is nodeless. In 
the itinerant picture, the Fermi surface nesting between 
the hole and the electron pockets facilitates the antifer- 
romagnetic fluctuations which favor the S2.2j,2-pairingi^. 
In real space, an intuitive picture of the Sa,2y2-wave pair- 
ing is just the next-nearest-neighbor (NNN) spin-singlet 
pairing with the s-wave symmetry^. Because the anion 
locations are above or below the centers of the iron-iron 
plaquettes, the NNN antiferromagnetic exchange J2 is at 
the same order of the nearest neighbor (NN) one Ji . The 
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NNN Sj.2j,2-wave pairing can be obtained from the decou- 
pling of the J2-term. On the other hand, various exper- 
imental results have shown signatures of nodal pairing 
structureai^^— . In the framework of the s-wave pairing, 
the nodal pairing can be achieved through the s,j.2_f_y2- 
pairingi^^— . The possibility for the Sa;2+j^2-wave pairing 
in iron-based superconductors has also been shown in the 
functional renormalization-group calculatioi*^^. 

Another important aspect of the iron-based supercon- 
ductors is the spontaneous anisotropy of both the lat- 
tice and electronic degrees of freedom, which reduces 
the 4-fold rotational symmetry to 2-fold. For example, 
LaOFeAs undergoes a structural orthorhombic distortion 
and a long-range spin-density wave (SDW) order at the 
wavevector (7r,0) or (0,7r)^. Similar phenomenon was 
also detected in NdFeAsO by using polarized and unpo- 
larized neutron-diffraction measurements-^. One popu- 
lar explanation of the nematicity is the coupling between 
lattice and the stripe-hke SDW oidev^^. The SDW or- 
dering has also been observed in the FeTe system but 
with a different ordering wavevector at (|-, 

Very recently, the experimental results of the FeSe 
superconductor reported by Song et al} indicate a 
pronounced nodal pairing structure in the scanning- 
tunneling-spectroscopy. Strong electronic anisotropy is 
observed through the quasi-particle interference of the 
tunneling spectra at much higher energy than the su- 
perconducting gap. The low energy tunneling spectra 
around the impurity and the vortex core also exhibit the 
anisotropy. The shapes of the vortex cores are signifi- 
cantly distorted along one lattice axis. 

The anisotropy may arise from the structural tran- 
sition from tetragonal to orthorhombic phase at 90K. 
However, the typical orthorhombic lattice distortions in 
iron superconductors are at the order of 0.0 12 A, which 
is about half a percent of the lattice constant and only 
lead to a tiny anisotropy in electronic structuresi. In 
contrast, the anisotropic vortex cores and impurity tun- 
neling spectra observed by Song et al. are clearly at 
the order of one. Therefore these anisotropics should be 
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mainly attributed to the electronic origin. The antifer- 
romagnetic long-range-order in such a system may be a 
possible reason for the anisotropy. For example, it has 
been theoretically investigated that the stripe-like SDW 
order can induce strong anisotropy in the quasi-particle 
interference of the STM tunnelling spectroscopy^^. How- 
ever, no evidence of magnetic ordering has been found in 
FeSe at ambient pressure^®, thus this anisotropy should 
not be directly related to the long-range magnetic order- 
ing. 

On the other hand, orbital ordering is another pos- 
sibility for nematicity in transition metal oxides. For 
example, orbital ordering serves as a possible mecha- 
nism for the nematic metamagnetic states observed in 
Sr3Ru207^^'^^, and its detection through quasi-particle 
interference has been investigated^Si^. Orbital ordering 
has also been suggested to lift the degeneracy between 
the dxz and dy^-orbitals to explain the anisotropy in iron- 
based superconductors. 

In this article, we will study the effect of orbital order- 
ing to the vortex tunneling spectra in the FeSe supercon- 
ductor. The rest of the paper is organized as follows. In 
Sec. ini a two-band model Hamiltonian and the relevant 
band parameters are introduced. The Bogoliubov-de 
Gennes mean-filed formalism is described in Sec. IIIII In 
Sec. IIVI we analyze the effects of orbital ordering to the 
tunneling spectra of mixed pairing of the NN 53,2 _|_j,2 -wave 
and the NNN S2;2j,2-wave in the homogeneous systems. In 
Sec. |Vl the effects of orbital ordering to the anisotropic 
vortex core tunneling spectra are investigated, which are 
in a good agreement with experiments. Discussions and 
conclusions are given in Sec. I VII 



II. MODEL HAMILTONIAN FOR THE BAND 
STRUCTURE 

For simplicity, we use the two-band model involving 
the dxz and dyz-orbital bands in a square lattice with 
each lattice site representing an iron atom, which was 
first proposed in Ref. (sTj . This is the minimal model de- 
scribing the iron-based superconductors, which can also 
support orbital ordering. The tight-binding band Hamil- 
tonian reads 



where denotes the creation operator for an elec- 
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a.<y,r — da ^ pda,a,r dcnotcs the particlc number oper- 



ator; ^ is the chemical potential; and t™ denote 
the longitudinal cr-bonding and transverse 7r-bonding be- 
tween nearest-neighboring (NN) sites, respectively; the 
three next-nearest-neighboring (NNN) hoppings can be 
expressed in terms of the NNN a and 7r-bondings i™" 

and i™", respectively as i^""" = ^(if"" + ^2"" = 

i(i|["" - i™"), and i^"" = 5(-i||"" + i™"). We de- 
pict the hopping schematic of the two-band tight-binding 
model Eq. ^ in Fig. [TJ 




FIG. 1: (Color online) The hopping schematic of the two-band 
tight-binding model Eq. U]) in a unit cell. Each black solid 
circle represents an Fe atom. The solid arrows denote NN lon- 
gitudinal (T-bonding (red) and transverse vr-bonding (green) 
hoppings, respectively. The NNN intra-orbital hoppings t""" 
are indicated by blue dot arrows along ±x 4- y directions. On 
the other hand, the NNN inter-orbital hoppings, tj"" and 
t"^"" , along the x + y and —x + y directions, respectively, are 
indicated by dash arrows labeled with cyan and pink. 



By introducing the spinor 'i'(fc) = [ilixz,a{k)^ i^yz^aik)]^ 
and performing the Fourier transformation, the Hamilto- 
nian in momentum space becomes 



(1) Ho^Y.'^l^^^)[Ha,{k)+e{k)-^l5ab}^bAk). (3) 



where 
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where a, b refer to the band index; the matrix kernel 
Hab{k) is written as 



2t"" cos kx + 2t™ cos ky 



4t3"" sin kx sin ky 
2<"" cos ky + 2i™ cos kx 



and e{k) = 2cosA:3, cosky[ff + i'^""). This two-orbital 
model has also been used to study impurity resonance 
states^^i^j vortex core states^i^, and quasiparticle scat- 
tering inference^. 



To fit the Fermi surface obtained from the LDA 
calculations^^, we use the parameter values below in the 
following discussions as 

if" = 0.8, tT = -1-4, tr" = 1-8, if"" = 0, (4) 

all which are in units of t = 1, which is roughly at the 
energy scale around lOOmeV. This set of hopping inte- 
grals shows a similar band structure, as the one Raghu 
et al. used. The band-width is about 14. In the following 
discussion, we use = 1.15 which corresponds to slightly 
hole doped regimes. The unfolded Brillouin zone (UBZ) 
embraces a hole surface around the F point [A; = (0,0)], 
and four hole pockets around the X point \k = (±7r, ±7r)] 
and four electron pockets around M point [k = (0, ±7r) 
or (±7r,0)^^. 

Our main purpose is to study the anisotropy effects in 
FeSe due to the orbital ordering. Orbital ordering has 
been proposed in iron-based superconductors in previous 
studies2§rd2.. Such an ordering may arise from the in- 
terplay between orbital, lattice, and magnetic degrees of 
freedom in iron superconductors. In this article, we are 
not interested in the microscopic mechanism of sponta- 
neous orbital ordering, but rather assume its existence to 
explain the vortex tunneling spectra observed in Ref. [l| . 
According to the experimental data^, strong anisotropy 
has already been observed at least at the energy scale of 
lOmeV, which is much larger than the pairing gap value 
around 2meV. Thus when studying superconductivity, we 
neglect the fluctuations of the orbital ordering, but treat 
it as an external anisotropy. For this purpose, we add an 
extra anisotropy term into the band structure Eq. ([l} as 

r.cr 

which makes the d^^^-orbital energy higher than that of 
dyz- For comparison, the Fermi surfaces without and 
with the anisotropy term Eq. are depicted in Fig. [5] 
(a) and (b), respectively. In Fig. [2](b), with the orbital 
term, the distortion of the electron and hole pockets in x 
and j/-directions of the Fermi surfaces appears such that 
the anisotropy is derived in the iron-based superconduc- 
tors. 



III. THE BOGOLIUBOV-DE GENNES 
FORMALISM 

In this section, we present the self-consistent 
Bogoliubov-de Gennes (BdG) formalism based on the 
band structure described in Sec. |lll In principle, for this 
multi-orbital system, the general pairing structure should 
contain a matrix structure involving both the intra and 
inter-orbital pairings. Here for simplicity, we only keep 
the intra-orbital pairing which is sufficient to describe the 
anisotropy observed in the experiment by Song et alX. 
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FIG. 2: Fermi surfaces in the UBZ for (a) = and (b) 
5e — 0.2. The horizontal and vertical axes denote k^c and 
ky, respectively. The parameter values used are if" = 0.8, 
tT = -1.4, tT" = 1.8, = 0, and fj, = 1.15. Anisotropic 
hole and electron pockets are shown in (b) at the F and M- 
points, respectively. 

The pairing interactions including the NN and NNN 
pairing are defined as 

H^nt = -| E E^a(^-',^")Aa(r,r'), 

(r,r') a 

- I E Y.^air,^)Aa{r,r^), (6) 

{{r,r')) a 

where a is the orbital index taking values of d^z and 
dyz', {r,'r') represents the NN-bonds and ((r, r*)) repre- 
sents the NNN-bonds; gi^2 denotes the pairing interac- 
tion strengths along the NN and NNN-bonds; Aa(r', r*) 
describes the spin singlet intra-orbital pairing operator 
across the bond defined as 

Aa(r,r') = da.,l.fda,-t,r' - da,t,rda.,l.f'- (7) 

where = r + 6. For the Sa;2_(_y2-pairing along the NN- 
bonds, S = ax{y), whereas for the Sa;2y2-pairing along 

the NNN-bonds, d = ±a{x + y) NNN-bonds, where a is 
Fe-Fe bond length, defined as the lattice constant. In 
the square lattice, these two pairings belong to the same 
symmetry class, thus they naturally coexist. After the 
mean-field decomposition, the Hamiltonian becomes 

Hmf = i?o-| E E^a(^^")Aa(r,r') 

(f.r') a. 

- f E E^a(^^")Aa(f,r') + /i.c., (8) 

{{r,r')) a 

where A*(r, r*) = (A]j(r, r*)) is the pairing order pa- 
rameter and (• • • } denotes the expectation value over the 
ground state. 

The mean-field BdG Hamiltonian Eq. can be diag- 
onalized through the transformation as 

fcaAr)\ f Ua,„(r) -<.„(r) \ /7a,n\ . . 
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The eigenvectors associated with En of the above BdG 
equations are (Ma,n(^); ''^a,n(^))"^ and the pairing order 
parameters can be further obtained self-consistently as 

Aa(r, r-) = J2{un{r)v:{r^) + Un{r^)vl{r)) tanh 

n 

(10) 

where b = l/fc^T. After the wave functions are obtained 
self-consistently, the local density of states (LDOS), 
which is proportional to the conductance (dl/dV) in the 
scanning tunneling microscopy, can be further measured 
by 

71, a 

+ \va,„{^\^L{E + En)Y (11) 

where L{x) is the Lorentzian as L{x) = 7/[7r(x^ + 7^)] 
and 7 means the energy-broadening parameters, usually 
set around 10~^. 

When we study the vortex lattice structure problem, 
the single particle Hamiltonian Eq. ([T]) is modified by 
the magnetic vector potential as 



(12) 



where $0 — hc/2e is the quantized flux; a, b denote or- 
bital indices and a denotes the spin index; ti,j.a,b repre- 
sents i™, i"" ti™3 depending on the corresponding 

bonds and orbitals. The vector potential A(f) is chosen 
as the Landau gauge by {Ax, Ay) = {0,Bx). 

Due to the magnetic translational symmetry, we can 
apply the magnetic periodic boundary conditions to form 
Abrikosov vortex lattices. Each magnetic unit cell car- 
ries a magnetic flux of 2<i>o, so that each magnetic unit 
cell contains two vortices. We choose the size of the mag- 
netic unit cell pa x qa with p — 2q, and the number of 
unit cells iV^; x Ny with Ny = 2Na:- The correspond- 
ing magnetic field is i? = 2^o/pqa^ — ^o/{qa)'^. In the 
Abrikosov vortex lattice, the translation vector is writ- 
ten as V = {Xpa,Yqa), where X = 0, • • • ,Nx — 1 and 
Y = 0, • • • , Ny — 1 are integers. The coordinate of an 
arbitrary lattice site can be expressed as .R = r + V, 
where r = {xa, ya) denotes the coordinate of the lattice 
site within a magnetic unit cell, i.e. 1 < a; < p and 

i<y<q- 

Under the magnetic periodic boundary conditions, the 
eigenvectors of the BdG equations Eq. ([9]) satisfy a pe- 
riodic structure written as 



Ua,n{r + pax) 
Va,nir + pax) 

Ua.n{^+ qay) \ 
Va,n{r + qay) J 



2-111— / ->\ 

e ''Ua,n(r) 

e-'"^i^a,„(r) 



,iK^ ( Ua,„ir) 
Va,n{r) 



Here = =^ and Ky = represent the magnetic 

Bloch wavevector on the x and y components'^'*. With 
the relation, we can simulate vortex lattices with sizes of 
{Nxpa) X (Nyqa) but reduce the computational effort by 
diagonalizing NxNy Hamiltonian matrices with dimen- 
sions of Apq rather than directly diagonalizing a ANxNypq 
Hamiltonian matrix. 



IV. THE NODAL V.S. NODELESS PAIRINGS 

In this section, we investigate the behavior of the su- 
perconducting gaps in the homogeneous system. Due to 
the translation symmetry, the pairing order parameters 
are spatially uniform, and we define lS.a{5) = Aa(r, f+5) 
where a = dxz, dyz- We start with the case in the absence 
of orbital ordering, i.e., 6e — 0. Due to the four- fold ro- 
tational symmetry, not all of the pairing order parame- 
ters are independent. For the NN bond pairing, we have 
^xzii) = \z{y) and Axziy) = \zix) due to the s- 
wave symmetry. As for the NNN-bonding, the similar 
analysis yields the following relations of Axz{S: + y) — 
Ayz[~x + y), and Axz{x - y) = Ay^{x + y). 

Due to the multi-orbital structure, generally speaking, 
the pairing order parameters have the matrix structure, 
thus the analysis of the pairing symmetry is slightly com- 
plicated. However, before the detailed calculation, we 
perform a simplified qualitative analysis by considering 
the trace of the pairing matrix, defined as 



1 



A{5)^-{Axz{5) + Ay^iS)) 



(14) 



for both of the NN and NNN-bonds. These quantities 
play the major role in determining the pairing symme- 
try. The angular form factor of the Fourier transform 
of the NN s-wave pairing is As^^^^^ (fc^, fc^) oc coskx + 
cos ky, and that of the NNN s-wave is As_^2 2{^x,ky) oc 
cos kx cos ky. Naturally As ^ ^ and A3,2j^2 have the same 
phase. Otherwise there would be a large energy cost cor- 
responding to the phase twist in a small length scale of 



lattice constant. The nodal lines of A, 



have intersec- 



(13) 



tions with the electron pockets, while the nodal lines of 
Ax2y2 have not intersections with Fermi surfaces. Gen- 
erally speaking, the NN and NNN s-wave pairings are 
mixed due to the same symmetry representation to the 
lattice group as 

Ag^ {kx, ky) — Ai(cos kx + cos ky) + A2 cos kx cos ky. 

(15) 

It is well-known that the gap function is nodal for the 
NN s-wave pairing; while it is nodeless for the NNN 
s-wave pairing. However, they can mix together. Re- 
cently this aspect was supported by a variational Monte- 
Carlo calculation,— where the authors discovered that 
the S2;2j,2-wave and Sa;2_|_j,2-wave states are energetically 
comparable. Our BdG calculations show that even for 
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the case of 52 7^ and gi = 0, the s-wave NN pairing is 
still induced by the 172 term, and vice versa for the s-wave 
NNN pairing showing that they can naturally coexist. 
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FIG. 3: DOS v.s. tunneling bias E for the mixed s-wave pair- 
ing states at the zero temperature. The red hollow squares 
with the parameter depict the nodeless pairing with the dom- 
inant NNN pairing; and the black solid line depicts the 
nodal pairing with the dominant NN pairing. The param- 
eter values are {gi = 0.8, 32 = 1-2) for the nodeless case and 
(gi — 1.2, (72 ~ 0.4) for the nodal case, respectively. 
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FIG. 4: (a) DOS v.s. tunneling voltage E for extended 
s±-wave without orbital anisotropy {Se — 0), and the or- 
bital anisotropy {Se = 0.2,0.3). Other parameter values are 
{gi — 1.2, g2 = 0.4). At & = the coherent peaks are around 
E = ±0.18 and upon increasing Se = the locations of the 
coherent peaks only slightly move towards zero energy, (b) 
The temperature dependence of the DOS v.s. E for the ex- 
tended s±-wave with orbital anisotropy Se = 0.2 and the same 
parameters of 511,2 as in (a). At T = 0, the coherent peaks are 
located at £ = ±0.165. 

For the coexistence of the NN and NNN s-wave pairing, 
the pairing gap function can be either nodal or nodeless. 
For the pure NNN s-wave pairing, the nodal lines of the 
A.,2,,2 are kx^y — which have no intersections with 
all of the hole and electron pockets, thus the pairing is 
nodeless. For the pure NN s-wave pairing, the nodal lines 
form a diamond box with four vertices at the M-points 
(±7r, 0) and (0, ±7r), which intersects both electron pock- 
ets, thus the pairing is nodal. When the NN and NNN 
s-wave pairings coexist, if the NN Sj.2+j,2-wave pairing 



is dominant, the nodal lines of the pairing function is 
illustrated in the Fig. S6 (a) of Ref. The origi- 

nal diamond nodal box is deformed by pushing the four 
vertices way from the M-points to the direction of the 
F-point. If the deformation is small, the deformed dia- 
mond still intersects with the electron pockets, and thus 
the pairing remains nodal. Upon increasing the strength 
of the NNN pairing, the deformation is enlarged and the 
intersections disappear. Thus the pairing is nodeless. 

Fig. [3] reveals the density of states (DOS) v.s. tunnel- 
ing voltage for the mixed s-wave pairing state. The red 
hollow squares are obtained using stronger NNN pairing 
strength [i.e. gi = 0.8 and g2 = 1-2 in Eq. providing 
gapful behavior which is similar to the pure Sx^y^-'^B.ve 
state. In this case, Axz{x) = Ayz{y) = 0.068, Axz{y) — 
Ayz{x) = -0.077 and Axz,yz{±i ± 2/) = 0.085 showing 
stronger NNN S2.2j,2-wave pairing than NN S2,2_|.j,2-wave 
pairing. On the other hand, the black solid line (using 
gi = 1.2 and g2 = 0.4) indicates a gapless ^-shape, sim- 
ilar to the pure Sa;2_|_j,2-wave state. In contrary to the 
former case, the NN pairing order parameters A^zix) — 
Ayziv) = 0.071, Axziy) = Ayz{x) = -0.059 are larger 
than NNN S3,2j,2-wave ones Axz,yz{^x±y) = 0.057. This 
reveals that the competition between the gapful and gap- 
less modes can be adjusted by tuning the ratios between 
NNN and NN pairing interactions. 

The main goal of our paper is the effect from or- 
bital ordering which is mimicked by Eq. ([5|). With the 
anisotropy, the pairing structure changes from Eq. (jl5p 
to the following 

As^{kxTky) — Ai (cos fcj; + A cos /cj,) 

+ A2 cos kx cos ky , (16) 

where A is determined by the anisotropy. The anisotropic 
nodal curve still intersects the electron pockets as illus- 
trated in the Fig. S6 (b) in Ref. 'd^ and thus the super- 
conducting gap functions remain nodal. The calculated 
DOS vs energy patterns are plotted in Fig. |l](a) with the 
parameter values specified in the figure caption. Upon fi- 
nite fc, the pairing order parameters become anisotropic. 
For example, at Se = 0.2, we have 

Axz{x) = 0.059, Ayziv) = 0.070, 
A:..(y) - -0.050, Ayzix) ^ -0.05A, 
Axzix + y) = Axz{~x + y) = 0.0503 
Ayzix + y) = Ayzi-x + y)= 0.0505, (17) 

and the gapless gap function and the V-shaped spectra 
remain in the moderate anisotropy from orbital ordering. 

We also check the temperature dependence of DOS 
with the orbital ordering as presented in Fig. 2] (b) with 
same parameter values of gi^2 as in Fig. 0] (a). With 
orbital anisotropy 6e = 0.2, the coherence peaks at zero 
temperature is approximately Ach ~ 0.165. At low finite 
temperatures, say T = 0.05 « 0.3Ach, the V-shape DOS 
is still discernible. However, upon increasing tempera- 
tures the y-shaped LDOS patterns smear and eventually 
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the coherent peaks disappear at T = 0.1 w 0.6Ach- In 
this case the system turns into the normal state. This 
feature is quahtatively consistent with the experimen- 
tal observation of the differential conductance spectra on 
FeSe in Ref. [1]. 



V. THE VORTEX STRUCTURE 

In this section, we will study the vortex tunneling spec- 
tra for the extended s-t-wave state. The size of the 
magnetic unit cell is chosen as pa x qa — 20a x 40a, 
which contains two vortices. The external magnetic field 
B — 2$o/W2^- The number of magnetic unit-cells shown 
below is using x Ny — 20 x 10, which is equivalent to 
the system size of 400a x 400a. The BdG equations are 
solved self-consistently with the tight-binding model Eq. 
(IT^ plus the mean-field interaction Eq. ©. The vortex 
configurations are investigated for both cases with and 
without orbital ordering in Sec. IV Al and Sec. IV Bl re- 
spectively. The interaction parameters are gi = 1.2 and 
g2 = 0.4, and the temperature is fixed at zero. 



A. Vortex structure in the absence of orbital 
ordering 

We start with the vortex configuration without or- 
bital ordering. The NN S3;2+j^2-wave and NNN Sx2y2- 
wave pairing order parameters in real space as defined as 
follows. We define the longitudinal and transverse NN 
s-wave pairings as 

A^^(f) = ^\^A,,{f,r + ax) + A,,{r,r~ ax) 

+ Ay^Xr,r + ay) + Ayz{f,f- ay)^ 

A^^{r) = ^i^Ay,ir,r + ax) + Ay,{r,f^ax) 

+ A^,{f,f+ay) + A^Af,f~ay)y (18) 
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FIG. 5: The vortex structure for the mixed NN s^2^y2-wave 
NNN 83.2^2 -wave pairing witliout orbital ordering, (a) and (b) 
depict the spatial distribution of the longitudinal and trans- 
verse NN pairings A^^(r^ and Ay^(r^, respectively, (c) de- 
scribes the NNN s-wave pairing order parameters A^^^(r). 
(d) The LDOS v.s. E at different locations from the vortex 
core [at r = (0, 10a)] to outside along the x-axis. The dis- 
tances of each site from the vortex core take the step of one 
lattice constant, (e) The spatial LDOS distribution at the 
energy of the vortex core resonance peak Ere ~ 0.062. 



For the NNN pairing related to site r, we define 

A™(f) = i E Aa{r,r + S). (19) 

a—xz,yz-^5—ixiy 

The pure NN 5^2 _|_j,2 -wave vortex states are recently in- 
vestigated to describe the competition between the super- 
conductivity and spin-density-wave (SDW) in the hole- 
doped materials, Bai„xKa;Fe2As2^. Similar calculation 
for the pure NNN 53,2 j^2 -wave with SDW has been used 
to study BaFei_a;Co2;As2'^^ and the FeAs stoichiometric 
compounds^. In our case, these two s-wave pairing order 
parameters mix together. 

The real space profiles of the longitudinal and trans- 
verse NN s-wave pairings A^^{r) and A^'^{r) are de- 
picted in Fig. [5] (a) and (b), respectively. The vortex 



cores are located at r = (0, ±10a) where the pairing or- 
der parameters are suppressed. Both of them exhibit the 
C4 symmetry. The NNN s-wave pairing order parameters 
are depicted in Fig. [S] (c) and the vortices are diamond- 
shaped with the C4 rotational symmetry. Note that the 
maximum magnitudes of the NNN S3;2j,2-wave pairing or- 
der parameters are smaller than those of NN longitudi- 
nal and transverse Sx2_^_y2-wave pairing order parameters, 
due to the stronger NN pairing strength {gi = 1.2 and 
52 = 0.4). The coherence length can be estimated as 
^ w 4a 5a from the spatial distributions of order pa- 
rameters. 

The relations of LDOS v.s. the tunneling energy E 
are presented in Fig. ^{d) at different locations from the 
vortex core to outside along a;-axis. The LDOS pattern 
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along y-axis is the same as Fig. [5](d) due to C4 symmetry. 
Note that there exist fine osciUations in the LDOS pat- 
tern. This fine oscillation structure may come from the 
Landau oscillation in which the oscillation period is re- 
lated to the external magnetic ffeld"^"* . At the vortex core 
[r = (0, 10a)], the coherence peaks at the bulk gap value 
of Ach disappear. Instead, a resonance peak appears 
at Ere ~ 0.062. Away from the vortex core, the reso- 
nance peak splits into the particle and hole branches with 
energies which symmetrically distribute with respect to 
Ef. As the distance increases, the peak intensities de- 
crease, and the energy separations between the particle 
and hole branches of peaks increase. As the distance 
reaches around 6a, i.e. beyond the coherence length ^, 
these peaks merge into the bulk coherence peaks. Fig. [S] 
(e) presents the spatial distribution of the LDOS at the 
vortex core resonance state energy Ere, which exhibits 
the 4-fold rotational symmetry. The vortex core state 
mainly distributes within one coherence length, thus it is 
closer to a bound state rather than a resonance state. 



B. The vortex structures with orbital ordering 
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FIG. 6: The spatial distributions of the s-wave pairing order 
parameters with the orbital anisotropy e = 0.2. (a) The lon- 
gitudinal NN pairing (r), (b) the transverse NN pairing 
Ay'^(r), and (c) the NNN s^2y2-weive pairing order parame- 
ters A^'^^ (r) . AH of them show anisotropy. 



In this subsection, we consider the effect of orbital or- 
dering to the vortex lattice states. The band and interac- 
tion parameters are the same as in Sec. IV Al except that 
we add the anisotropy term of Eq. ([5]) with 5e — 0.2. 
Such a term breaks the degeneracy between the d^z and 
dyx-orbitals, and reduces the C4 symmetry down to C2. 

Fig. M (a), (b) and (c) depict the spatial distribu- 
tions of the NN Sa;2+y2 -pairing order parameters A;^^ 
and Ay^, and the NNN 53.2^^2 -pairing order parameters 
in a magnetic unit cell, respectively. All of them clearly 
exhibit the breaking of the 4-fold symmetry down to the 
2-fold one. For the dominant NN S2,2^j^2-pairings, the 



coherence lengths along the x and y-directions are no 
longer the same, which can be estimated as ~ 3a and 
^j, w 7a, respectively. From Fermi surface Fig. [2](b), in 
the presence of orbital ordering, the electron pocket in 
the y-direction shrinks, which implies that Cooper pair- 
ing superfluid stiffness is weaker along the y-direction 
than the x-direction. This picture agrees with the larger 
value of £,y exhibiting in Figs. [SI 
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FIG. 7: The LDOS v.s. tunneling energy E at different sites 
from vortex core to outside along (a) the a;-axis (short axis) 
and (b) the y-axis (long axis). The distances of each site from 
the vortex core take the step of one lattice constant, (c) The 
spatial LDOS distributions at the vortex core resonant energy 

Eer = 0.038 . 

Next we turn to study the LDOS patterns for the ex- 
tended s-|--wave with orbital ordering. In comparison 
with those without orbital ordering depicted in Fig. [S] 
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(c), the LDOS patterns in Fig. [7] (a) and (b) exhibit 
significant anisotropy. At the vortex core center, the 
resonance peak sphts into two pieces inside gap. This 
feature is similar to the previous studies on the s^iyi- 
wave pairing with SDW in iron-based superconductor's^ 
and the d-wave pairing with antifcrroniagnetic ordering 
in cuprates^. Away from the vortex center, however, 
the peaks of the particle and hole branches along short 
(x) and long {y) axes behave differently. The separation 
between the particle and hole peaks disperses along the 
short axis much quicker than that along the long axis. 
The peak intensities along the short axis are stronger 
than those along the long axis. These features are in 
good agreement with recent experiment observations^. 
The pronounced anisotropy also exhibits in the spatial 
variation of the LDOS at the energy E^r — 0.038 where 
the resonance peak is located, presented in Fig. [7] (c). 
Moreover, with negative which loads particles in dxz 
in prior to dyz, all above pictures will have a 7r/2-rotation. 
Therefore, the experimentally observed anisotropy agrees 
with the picture of orbital ordering. 

VI. DISCUSSIONS AND CONCLUSIONS 

We have studied a minimal two-orbital model with or- 
bital ordering to interpret the recent STM observations 
in Ref. [l| including the nodal superconductivity and the 
anisotropy vortex structure in FeSe superconductors. In 
considering the absence of magnetic long-range-order in 
FeSe at ambient pressure, orbital ordering provides a nat- 



ural formalism for anisotropy. The NN s^+y^i-^SMe and 
the NNN s^2y2-wave are considered, which generally are 
mixed. When the NNN pairing dominates, nodal pair- 
ing still exists even in the appearance of orbital ordering. 
We further performed the BdG calculation for the vor- 
tex tunneling spectra in the presence of orbital ordering, 
which breaks C4 rotational symmetry down to C2 and 
explicitly induces anisotropic vortex structures. 

The microscopic mechanism of the origin of this spon- 
taneous anisotropy remains an open question. It might 
be related to the strong antiferromagnetic fluctuations. 
As shown in Ref. [2Tll22ll50j . before the onset of the an- 
tiferromagnetic long rang order, the spin nematic order 
with a Z2 symmetry breaking occurs. This nematic order 
corresponds to anisotropic spin-spin correlation along a 
and 6-axis, which can in turn induce orbital ordering. 

VII. ACKNOWLEDGEMENT 

C. W. thanks Z. Y. Lu, and F. J. Ma for help- 
ful discussion. H. H. H. and C. W. are partly sup- 
ported by the NBRPC (973 program) 2011CBA00300 
(2011CBA00302), NSF-DMR-1105945 and Sloan Re- 
search Foundation. C. L. S., X. C, X. C. M. and Q. 
K. X. are supported by National Science Foundation and 
Ministry of Science and Technology of China. 

Note added Near the completion of this manuscript, 
we learned the article by Chowdhury et al.^'^ which stud- 
ied the anisotropic vortex tunneling spectra in Ref. [l| 
through the Ginzburg-Landau formalism. 



^ Can-li Song, Yi-lin Wang, Peng Cheng, Ye-ping Jiang, Wei 
Li, Tong Zhang, Zhi Li, Ke He, Lili Wang, Jin-feng Jia, 
Hsiang-Hsuan Hung, Congjun Wu, Xucun Ma, Xi Chen, 
and Qi-kun Xue, Science 332, 1410 (2011). 

^ Y. Kamihara, T. Watanabe, M. Hirano and H. Nosono, J. 
Am. Chem. Soc. 130, 3296 (2008). 

^ C. de la Cruz et al, Nature (London) 453, 899 (2008). 

* T. Hanaguri, S. Nittaka, K. Kuroki and H. Takagi, Sceince, 
328 (2010). 

^ X. Zhang, Y. S. Oh, Y. Liu, L. Yan, K. H. Kim, R. L. 
Greene, and I. Takeuchi, Phys. Rev. Lett. 102, 147002 
(2009). 

H. Ding et al, Europhys. Lett. 83, 47001 (2008). 

K. Seo, B. A. Bernevig and J. Hu, Phys. Rev. Lett. 101, 

206404 (2008). 

^ M. Daghofer et al, Phys. Rev. Lett. 101, 237004 (2008); 

A. Moreo, M. Daghofer, J. A. Riera, and E. Dagotto, Phys. 

Rev. B 79, 134502 (2009). 
^ F. Wang, H. Zhai, Y. Ran and A. Vishwanath, D. H. Lee, 

Phys. Rev. Lett., 102, 047005 (2009). 

V. Cvetkovic and Z. Tesanovic, Europhysics Lett. 85, 

37002 (2009). 

V. Cvetkovic and Z. Tesanovic, Phys. Rev. B 80, 024512 
(2009). 

M. R. Norman, Physics 1, 21 (2008). 



J. D. Fletcher et al, Phys. Rev. Lett. 102, 147001 (2009). 
C. W. Hicks et al, Phys. Rev. Lett. 103, 127003 (2009). 
^5 B. Zeng, G. Mu, H. Q. Luo, T. Xiang, H. Yang, L. Shan, 
C. Ren, 1. 1. Mazin, P. C. Dai and H.-H. Wen, Nat. Comm. 
1, 112 (2010). 

^® I.I. Mazin, D.J. Singh, M.D. Johannes and M.H. Du, Phys. 

Rev. Lett. 101, 057003 (2008). 
" Z.-J. Yao, J.-X. Li and Z. D. Wang, New. J. Phys. 11, 

025009 (2009). 

F. Wang, H. Zhai and D. -H. Lee, Phys. Rev. B 81, 184512 
(2010). 

Fa Wang, Hui Zhai, Ying Ran, Ashvin Vishwanath and 
Dung-Hai Lee, Phys. Rev. Lett. 102, 047005 (2009). 
^° Y. Chen et al, Phys. Rev. B 78, 064515 (2008). 

Chen Fang, Hong Yao, Wei-Feng Tsai, Jiang Ping Hu and 
Steven A. Kivelson, Phys. Rev. B 77, 224509 (2008). 
Cenke Xu, Markus MuUer and Subir Sachdev, Phys. Rev. 
B 78, 020501(R) (2008). 

Shiliang Li, Clarina de la Cruz, Q. Huang, Y. Chen, J. W. 
Lynn, Jiangping Hu, Yi-Lin Huang, Fong-chi Hsu, Kuo- 
Wei Yeh, Maw-Kuen Wu and Pengcheng Dai, Phys. Rev. 
B 79, 054503 (2009). 
^'^ Wei Bao, Y. Qiu, Q. Huang, M.A. Green, P. Zajdel, M.R. 
Fitzsimmons, M. Zhernenkov, S. Chang, M. Fang, B. Qian, 
E.K. Vehstedt, J. Yang, H.M. Pham, L. Spina and Z.Q. 



9 



Mao, Phys. Rev. Lett. 102, 247001 (2009). 

J. Knolle, I. Ercmin, A. Akbari, and R. Moessner, Phys. 

Rev. Lett. 104, 257001 (2010) 
^'^ S. Medvedev et al, Nat. Mater. 8, 630 (2009). 

W. C. Lee and C. Wu, Phys. Rev. B 80, 104438 (2009). 
^® S. Raghu, A. Paramekanti, E-.A. Kim, R. A. Borzi, S. A. 

Grigera, A. P. Mackenzie and S. A. Kivelson, Phys. Rev. 

B 79, 214402 (2009). 

W. C. Lee, D. P. Arovas and C. Wu, Phys. Rev. B 81, 
184403 (2010) 

^° W. C. Lee and C. Wu, Phys. Rev. Lett. 103, 176101 (2009). 
3^ S. Raghu, X. -L. Qi, C. -X. Liu, D. J. Scalapino and S. -C. 

Zhang, Phys. Rev. B 77, 220503(R) (2008). 
32 W. -F. Tsai, Y. -Y. Zhang, C. Fang and J. Hu, Phys. Rev. 

B 80, 064513 (2009). 

T. Zhou, X. Hu, J. -X. Zhu and C. S. Ting, e-print: 
0904.4273. 

^'^ X. Hu, C. S. Ting and J.-X. Zhu, Phys. Rev. B 80, 014523 
(2009). 

H. -M. Jiang, J.- X. Li and Z. D. Wang, Phys. Rev. B 80, 

134505 (2009). 

Y. -Y. Zhang et al, Phys. Rev. B 80, 094528 (2009). 

S. Lebegue, Phys. Rev. B 75, 035110 (2007); D. H. Lu et 

al, Science 455, 81 (2008). 
3« Frank Kriiger et al, Phys. Rev. B 79, 054504 (2009). 
39 C. -C. Lee, W. -G. Yin and W. Ku, Phys. Rev. Lett. 103, 



267001 (2009). 

*° Wcichcng Lv, Jianshcng Wu, and Philip Phillips, Phys. 
Rev. B 80, 224506 (2009). 

C. -C. Chen, J. Maciejko, A. P. Sorini, B. Moritz, R. R. P. 
Singh and T. P. Devereaux, Phys. Rev. B 82, 100504(R) 

(2010) . 

Weicheng Lv, Frank Kriiger and Philip Phillips, Phys. Rev. 
B 82, 045125 (2010). 
"3 Wcichcng Lv and Philip Phillips, Phys. Rev. B 84, 174512 

(2011) . 

'^'^ Qiang Han, J. Phys.: Condens. Matter 22, 035702 (2010). 
F. Yang, H. Zhai, F. Wang and D. -H. Lee, Phys. Rev. B 
83, 134502 (2011). 

*® Can-li Song, Yi-lin Wang, Peng Cheng, Ye-ping Jiang, Wei 
Li, Tong Zhang, Zhi Li, Kc He, Lili Wang, Jin-fcng Jia, 
Hsiang-Hsuan Hung, Congjun Wu, Xucun Ma, Xi Chen, 
and Qi-kun Xue, Supporting Online Material in Science 
332, 1410 (2011). 

Y. Gao et al, Phys. Rev. Lett. 106, 027004 (2011). 
J.-X. Zhu and C. S. Ting, Phys. Rev. Lett. 87, 147002 
(2001). 

D. Chowdhury, E. Berg, and S. Sachdev, Phys. Rev. B 84, 
205113 (2011). 

^° R. M. Fernandes, A. V. Chubukov, J. Knolle, and I. 
Eremin, and J. Schmalian, Phys. Rev. B 85, 024534 (2012) 



